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We study second-order cosmological perturbations in scalar-tensor models of dark energy that 
f*^ , satisfy local gravity constraints, including f(R) gravity. We derive equations for matter fluctuations 

^sj ■ under a sub-horizon approximation and clarify conditions under which first-order perturbations in 

the scalar field can be neglected relative to second-order matter and velocity perturbations. We 
also compute the skewness of the matter density distribution and find that the difference from the 
ACDM model is only less than a few percent even if the growth rate of first-order perturbations is 
significantly different from that in the ACDM model. This shows that the skewness provides a model- 
independent test for the picture of gravitational instability from Gaussian initial perturbations 
including scalar-tensor modified gravity models. 
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I. INTRODUCTION 



The constantly accumulating observational data [l| continue to confirm that the Universe has entered the phase 

of an accelerated expansion after the matter-dominated epoch. To reveal the origin of dark energy (DE) responsible 

for this late-time acceleration is one of the most serious stumbling block in modern cosmology [2], |3( • The first step 

r~~^ toward understanding the nature of DE is to find a signature whether it originates from some modification of gravity 

t-H , or it comes from some exotic matter with negative pressure. If gravity is modified from Einstein's General Relativity, 

^D ' this leaves a number of interesting experimental and observational signatures that can be tested. Especially local 

gravity experiments generally place tight bounds for the parameter space of modified gravity models. 

[ — ■ So far many modifieclgravity DE models have been proposed- ranging from f(R) gravity [4[ (R is a Ricci scalar), 

scalar-tensor theory [5|, [6[ to braneworld scenarios [7J. The f(R) gravity is presumably the simplest generalization to 

the A-Cold Dark Matter (ACDM) model (f(R) = R— A). Nevertheless it is generally not easy to construct viable 

f(R) models that satisfy all stability, experimental and observational constraints while at the same time showing 

J> , appreciable deviations from the ACDM model. In order to avoid that a scalar degree of freedom (scalaron) as well 

as a graviton becomes ghosts or tachyons we require the conditions J,rr > and f.R > [8J- These conditions are 

m | also needed for the stability of density perturbations [9| . For the existence of a matter-dominated epoch followed by a 

late-time acceleration, the models need to be close to the ACDM model (m = Rf,RR,/f,R ~ +0) in the region R^ Ro 

(Ro is the present cosmological Ricci scalar) [To| ■ Moreover the mass of the scalaron field in the region R 3> Rq is 

sufficiently heavy for the compatibility with local gravity experiments [ll|, LLfi [Dl, H| • Finally, for the presence of 

a stable de Sitter fixed point at r = —Rf t R,/f = —2, we require that < m(r — —2) < 1 [10, Ha] . The models 

proposed by Hu and Sawicki [161 ] and Starobinsky 8] satisfy all these requirements. They take the asymptotic form, 

f(R) ~ R — uRrfl — (R/ R r ) ~ 2n ] (/i > 0, R c > 0, n > 0), in the region i? > R c (R c is roughly the same order as Ro). 

See Refs. [H Dl Hfl El for other viable -^ models - 

The main reason why viable f(R) models are so restrictive is that the strength of a coupling Q between dark energy 

and non-relativistic matter (such as dark matter) is large in the Einstein frame {Q = — l/v6) [22J. In the region of 

high-density where local gravity experiments are carried out, the scalaron field <f> needs to be almost frozen [IJ, [l2| 

with a large mass through a chameleon mechanism [23J] to avoid that the field mediates a long ranged fifth force. 



'Electronic address: 
t Electronic address: 



tatckawa@cpd.kogakuin.ac.jp 
shinji@rs.kagu.tus.ac.jp 



Cosmologically this means that the field does not approach a kinematically driven matter-dominated era ( "0MDE" 
|24 |) in which the evolution of scale factor is non-standard (a oc t 1 / 2 [22]). The deviation from the ACDM model 
becomes important as the mass of the scalaron gets smaller so that the field begins to evolve slowly along its potential. 
In other words the effect of modified gravity manifests itself from the late-time matter era to the accelerated epoch 
[a Ell • This leaves a number of interesting observational signatures for the equation of state of DE [18|, [2(| , matter 
power spectra p, M, [20] and convergence spectra in weak lensing [23, [26| • 

One can generalize the analysis in f(R) gravity to the theories that have arbitrary constant couplings Q [27|. In 
fact this is equivalent to Brans-Dicke theory [28j] with a scalar-field potential V (</>). By designing the potential so 
that the field mass is sufficiently heavy in the region of high density, it is possible to satisfy both local gravity and 
cosmological constraints even when \Q\ is of the order of unity [23]. The representative potential of this type is g iven 
by V{<t>) = V Q \1 - C(l - e - 2W ) p ] (V > 0, C > 0, < p < 1), which covers the f(R) models of Hu and Sawicki pjf and 
Starobinsky [8| as special cases. Especially when \Q\ is of the order of unity, these models lead to the large growth 

of matter density perturbations (S oc tW 25+4S Q - 1 )/ 6 ) at a late epoch of the matter era compared to the standard 
growth (6 oc t 2 ' 3 ) at an early epoch. This gives rise to a significant change of the spectral index of the matter power 
spectrum relative to that in the ACDM model [8|, [20] . Moreover it was recently shown that the convergence power 
spectrum in weak lensing observations is subject to a large modification by the non-standard evolution of matter 
perturbations [251 ]. 

In this paper we shall study another test of modified gravity DE models mentioned above by evaluating a normalized 
skewness, S3 = (S 3 )/(S 2 ) 2 , of matter perturbations. The skewness provides a good test for the picture of gravitational 
instability from Gaussian initial conditions [29j | . If large-scale structure grows via gravitational instability from 
Gaussian initial perturbations, the skewness in a Universe dominated by a pressureless matter is known to be S3 = 34/7 
in General Relativity [30 ]. Even when cosmological constant is present at late times, the skewness depends weakly 
on the expansion history of the Universe (less than a few percent) [33, [32|, 1331 - This situation hardly changes in 
open/closed universes [34[| and Dvali-Gabadadze-Poratti braneworld models [35[. One can see some difference for 
the models that are significantly different from Einstein gravity-such as Cardassian cosmologies [35l l36j . modified 
gravity models that respect Birkhoff's law [33]. In the context of dark energy coupled with dark matter, it was 
shown in Ref. [3a ] that the skewness can be a probe of the violation of equivalence principle between dark matter and 
(uncoupled) baryons. 

In Brans-Dicke theory with cosmological constant A the skewness has been calculated in Ref. [39] under the condition 
that the Brans-Dicke field is massless. In this case the evolution of scale factor during the matter-dominated epoch 
is given by a(t) oc £( 2 "bd+2)/(3w B d+4) jjj^ where t is a cosmic time and ojbd is a Brans-Dicke parameter. If the 
field is massless, the Brans-Dicke parameter is constrained to be wbd > 40000 [40] from solar-system experiments. 
This shows that the evolution of the scale factor in the matter era is very close to the standard one: alt) oc t 2 / 3, . 
We note that an effective gravitational "constant" that appears as a coefficient of matter density perturbations is 
also subject to change in Brans-Dicke theory. However it was found that the skewness in such a case is given by 
S3 = (34cjbd + 56)/(7wbd + 12) _39J] during the matter era, which is very close to the standard one (S3 = 34/7) under 
the condition wbd > 40000. 

The f(R) gravity corresponds to theory with the Brans-Dicke parameter wbd = [41]. Even in this situation, if 
the scalaron field has a potential whose mass is sufficiently large in the region of high density, the f(R) models can 
pass local gravity constraints as in the models proposed in Refs. [8|, [l6[. In such cases, compared to Brans-Dicke 
theory with a massless field, it is expected that the skewness may show significant deviations from that in General 
Relativity. Since the evolution of scale factor and matter perturbations is different from that in the massless case, we 
can not employ the result of the skewness presented above. 

In this paper we study second-order perturbations and the skewness for Brans-Dicke theory in the presence of a 
potential V(4>). This is equivalent to the scalar-field action given in Eq. ([2]) by identifying the coupling Q with the 
Brans-Dicke parameter ojbd via the relation 1/(2Q 2 ) = 3 + 2ojbd- In the massless case the solar-system constraint, 
^bd > 40000, gives the bound \Q\ < 10~ 3 , but it is difficult to find some deviations from General Relativity in such 
a situation. Our interest is the case in which the coupling Q is of the order of 0.1 < |Q| < 1 with the field potential 
that has a sufficiently large mass in the high-density region. This analysis includes viable f(R) models [8|, [l6| recently 
proposed in the literature. We would like to investigate how much extent the skewness differs from that in the ACDM 
model. We also derive conditions under which the contribution coming from first-order field perturbations can be 
neglected relative to second-order matter and velocity perturbations by starting from fully relativistic second-order 
perturbation equations. 



II. MODIFIED GRAVITY MODELS 



The action for Brans-Dicke theory [28[ in the presence of a potential V is given by 



S 



j4 / 

d Xy/-g 



1 „ W B D 



(Vx) 2 - V(x) 



+ Smig/tut^Fm) ■ 



(1) 



where x is a scalar field coupled to a Ricci scalar R, wbd is a so-called Brans-Dicke parameter and S m is a matter 
action that depends on the metric p M „ and matter fields ^ m . We shall use the unit 8irG = 1, but we restore the bare 
gravitational constant G when it is required. 



The action {1} is equivalent to the following scalar-tensor action with the correspondence x 



-iQ<i>. 



S = 



d i x^ 



where 



iF(0)i?-^(0)(V0) 2 -V(0) 



+ S m (g^,^m) , 



(2) 



f(0)=e- 2C W, u(4>) = (1 - 6Q 2 )F(0) . (3) 

As we already mentioned, the constant Q is related with wbd via the relation 1/(2Q 2 ) = 3 + 2wbd- In the limit Q — > 
(i.e., ojbd — > oo), the action @ reduces to the one for a minimally coupled scalar field (f> with a potential V(<^)- The 
f(R) gravity corresponds to the coupling Q = — lVo, i.e., wbd = 0. 

In the absence of the potential V((f>) the coupling Q is constrained to be \Q\ < 10~ 3 from solar-system tests. We 
are interested in the case where the presence of the potential can make the models be consistent with local gravity 
constraints (LGC) even for \Q\ — 0(1). The representative potential of this type is given by [27| 

V{<j>) = V [1 - C(l - e- 2W f] (Vo>Q,C>0,0< P <l), (4) 

where Vq is of the order of the present cosmological Ricci scalar Rq in order to be responsible for the acceleration of 
the Universe today. Note that the f(R) models proposed by Hu and Sawicki [16[ and Starobinsky .8] take the form 
f(R) = R— f-iR c [l - (R/R c )~ 2n \ (ji > 0, R c > 0, n > 0) in the region R > R c . These f(R) models are covered in the 
action © with (jj) by identifying the field potential to be V = (RF - /)/2 with F = df/dR = e 2 ^/^. 

The background cosmological dynamics and LGC for the potential ^ have been discussed in details in Ref. [27J . 
We review how the matter-dominated era is followed by the stage of a late-time acceleration. This is important 
when we discuss the evolution of matter density perturbations in Sec. IIV1 As a matter source we take into account a 
non-relativistic matter with an energy density p m . In the flat Friedmann-Lemaitre-Robertson- Walker (FLRW) metric 
with scale factor a(t), where t is cosmic time, the evolution equations for the action © are 



3FH 2 = -u(j} z +V- 3HF + p m , 



2FH 



-ojc/) 2 



F + HF 

F 



p m + 3Hp m = , 



V t 4--F,4,R = 0, 



(5) 
(6) 

(7) 

(8) 



where H = a/a is the Hubble parameter and a dot represents a derivative with respect to t. Note that the Ricci 
scalar is given by R — 6(2H 2 + H). 

We introduce the following dimcnsionless quantities: 



Xi 




and 



^ Z'rn — 



3FH 2 



= 1 - (1 - QQ 2 )x\ -x\- 2V6Q Xl , 



where we used Eq. ([5]). We then obtain 
dxi V6 



d^V 
dx 2 
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where N = In (a) and A = — V^/V. The effective equation of state is defined by 

2 H 



Wcff 



-1 



3H 2 



where 



H 
IP 



1-6Q 2 



3 + ix\ - 3xj - QQ 2 x\ + 2V6Qa;i + ZQ(\x 2 2 - AQ) 



(13) 



(14) 



When A is a constant (i.e., V(4>) = Voe~ x ^), the fixed points of the system can be derived by setting dxi/dN = 
dx2/dN — 0. Even if A changes with time, as it is the case for the potential (|4]), the fixed points can be regarded as 
instantaneous ones. The following points can play the role of the matter-dominated epoch: 



(Ml) (f> matter-dominated era 



(xi,x 2 ) 



V6Q 



V 3(2Q 2 -1) : 
• (M2) "Instantaneous" scaling solution 

(xi,x 2 ) 



I l rn — 



3-2Q 2 
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(15) 
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2Q 
A 



(16) 



In order to realize the matter era (fi m ~ 1 and w e s — 0) by the point (Ml), we require the condition Q 2 < 1. This 
point was used in the coupled quintessence scenario [24| (in the Einstein frame) where the coupling is constrained to 
be |Q| ;$ 0.1 from Cosmic Microwave Background anisotropies. In f(R) gravity (Q = —1/V&) we have fl m = 2 and 
w c ff = 1/3 (i.e., a oc t 1 ^ 2 [22J]), which means that the point (Ml) can not be responsible for the matter era for |Q| of 
the order of unity. 

The matter era can be realized by the point (M2) for |A| ^ \Q\ = 0(1). The parameter A for the potential ([4]) is 
given by 



A: 



2CpQe- 2 Q' t >(l - , 



-2Q<j,\p-l 



1 - C(l - e- 2 Q*)P 



(17) 



which is much larger than 1 for \Q4>\ ^C 1 (provided that C and p are not very much smaller than 1). Since R ~ p m /F 
during the deep matter-dominated epoch, the field <j) is stuck at the instantaneous minima characterized by the 
condition V^(4> m ) + Qp m ~ [see Eq. (J7J]. For the potential J3]) this translates into 



1 ( 2V QP C 

2Q V Pm 



(18) 



which means that |Q0 TO | <§C 1 and hence |A| ~^> 1 during the deep matter era {p m ~^> Vq). When \Q\ = O(l), the 
matter era is realized by the point (M2) instead of (Ml). 

For the dynamical system given by Eqs. (fTTjl and (TT2|) there exist the following fixed points that lead to the late-time 
acceleration: 



• (Al) Scalar-field dominated point 
v / 6(4Q - A) 



(xi,x 2 ) 



6(4Q 2 -QA-1) : 



6 - A 2 + 8QA - 16Q 2 



6(4Q 2 -QA-1) 2 



1/2N 



a 



w cS 



20Q 2 - 9QX - 3 + A 2 
3(4Q 2 - QA - 1) 



.(19) 



(A2) de Sitter point (present for A = 4Q) 



(xi,x 2 ) = (0, 1) , Q m = 0, w e s = -l. 



(20) 



The de Sitter point (A2) appears only in the presence of the coupling Q (characterized by the condition V^ + QFR = 
in Eq. ([7]), i.e., A = 4Q). This can be regarded as the special case of the accelerated point (Al). For the potential (g]) 
the parameter |A| is much larger than \Q\ during the matter era, but it gradually becomes the same order as \Q\ as the 
system enters the accelerated epoch. It was shown in Ref. [27J| that the de Sitter point (A2) is stable for d\/d<j) < 0. 
As long as |A| continues to decrease with the growth of \<fi\, the solutions are finally trapped at the stable de Sitter 
point (A2). If the stability condition, d\/dcf> < 0, is not satisfied, the solutions approach another accelerated point 

In the following we are mainly interested in the case where the "instantaneous" matter point (M2) is followed by 
the de Sitter point (A2). During most stages of cosmic expansion history the field <p is trapped at instantaneous 
minima of an effective potential induced by the matter coupling. This means that the condition, 4> 2 -C H 2 , is well 
satisfied. 

The mass squared of the field cf> for the potential (j4|) is given by 

M 2 = V M = W CpQ 2 e- 2Q < t '(l - e -2Q^p-2 (1 _ pe -2Q^ } (2 1) 

which is much larger than i?o (~ Vo) in the reg ion R ^> i?o- In this situation it is possible to satisfy local gravity 
constraints in the region of high density [ll|, [l2|, LLI| through a chameleon mechanism [23[ . Since the field is massive 
inside a spherical symmetric body with radius r c , only the surface part of its mass distribution contributes to the 
field profile outside the body. The effective coupling Q c g between the field and the pressureless matter is suppressed 
by a thin-shell parameter Ar c /r c relative to the bare coupling Q. For the potential ((3]) it was shown in Ref. [12| 
that constraints coming from solar system tests as well as the violation of equivalence principle give the bounds 
p > 1 — 5/(9.6 — log 10 |<3|) and p > 1 — 5/(13.8 — log 10 |Q|), respectively. In f(R) gravity these constraints correspond 
to p > 0.50 and p > 0.65, respectively. 

Substituting the field value ([18]) for Eq. (|21[) . we find that the mass squared during the matter era is given by 



/ o2-p \ i- p / tt2\ i- p 

M ^(kc) (1 - p)02 (^) " 2 - (22) 

where we used 3H 2 ~ p rn . We then find that the inequality, M 2 ^> H 2 , holds for the values of p,C,Q not very 
much smaller than unity. In the next section we shall use this property when we derive the equation for matter 
perturbations approximately. 

III. SECOND-ORDER COSMOLOGICAL PERTURBATIONS 

In this section we consider second-order cosmological perturbations for the action ^ and derive the equation for 
matter perturbations approximately. 

A. Perturbation equations 

Let us start with a perturbed metric including scalar metric perturbations a, (3, ip and 7 about the flat FLRW 
background [42| |: 

d,s 2 = -(1 + 2a)dt 2 - 2af3. l dtdx i + a(t) 2 [(1 + 2</?)% + 2 7Ki ] da;W . (23) 

At the second-order the scalar variables are written as 

a = a«+a( 2 ), (3 = /3« + /?< 2 ) , p = tpM + <p™ , 7 = 7 « + 7 ^ , (24) 

where the subscripts represent the orders of perturbations. We introduce the following quantities: 

X = a(/3 + a 7 ), K = SK, (25) 

where SK is the perturbation of an extrinsic curvature K. 

We decompose the scalar field <f> and the quantity F into background and perturbed parts: 

4>=4> (t) + ^(t,x), F = F (t) + 6F(t,x), (26) 



where 8<f) and SF depend on t and a position vector x. In the following we omit the subscript "0" from background 
quantities. The components of the energy momentum tensor of a pressureless matter can be decomposed as 



Tq = -{p m + Sp m ) , Tf = -p m V,; 



(27) 



where v is a rotational-free velocity potential. At the second-order, the perturbed quantities can be explicitly written 
as 



= 5ft^+5ft 2 K SF = SFW+SFW, 8p m = 8p$ + 5pW, v = vV>+v<® 



(28) 



The perturbation equations for the action @, up to the second-order, have been derived in Ref. [43| (see also 
Refs. [44J]). They are given by 
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where 
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The equations for the perturbations 5<fi and SF are 



2V$ - FjR 

2co 



(24> + 3H^+^^ 



(35) 



(36) 
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ujft + 2F + 3HF) a + \f8R = N SF , 



(37) 



(38) 



where Ns,/, and Nsf are second-order terms whose explicit expressions are given in Ref. [43J | . 

At the first-order the quantity, Sp^ 1 ' = Spin — pm® 1 ' j is known to be gauge-invariant [4J]. In order to construct 
gauge-invariant variables at the second-order, we introduce the following quantities 

5p = 5p m - f>m,av + <5p (<z) , (39) 

v x = v - X - X + v® , (40) 

where 5p^ q > and v x are quadratic combinations of first-order terms. By defining 6 m = 5p m / p m , it was shown in 
Ref. [43| that the following quantity is gauge- invariant at the second-order: 

8 = 8 m + 3aHv - —av + ~p m Ha?v 2 - v>*v % - 3p m aH A _1 V* ( — v A , (41) 

Pm 2 \Pm ' J 

where 6p v = 8p m — p m o-v. Note that the quantity v x can be also made gauge- invariant [43j. 

B. Approximate second-order equations 

If we take the temporal comoving gauge (v = 0), we have S = 5 m and q : i — [see Eqs. (j2"?| and (|4"Tj) ]. Taking 7 = 
for the spatial gauge condition, it follows that (3 = x/ a from Eq. P5|) . From Eq. (I34[) we obtain 

a = -^f3^\ (42) 



which means that a is a second-order quantity. Up to the second-order, Eqs. (|31[) and (|33p are written as 



A _3 

H? X+ 2~F 
1 



k= — - x + —( u< l ) 8<l> + SF-H6F-Fa), (43) 



6-n = n6--5 l f3< 1 . (44) 

a : 

In order to evaluate the terms on the r.h.s. of Eq. (|4"4"1) , it is sufficient to consider Eqs. (|4"0|) and (|4"3")) at the first 
order with the gauge v = 0. We then have x = ~~ atl x' ^ = ~~ w x anc ^ 

ft (i) = _^x + — (u><j)d<f> + 5F-HSF), (45) 

a 2-F 



where we omitted the order of the subscript from the r.h.s. of these equations. Hence Eq. (J44j) can be read as 

5-n=-X7 ■ (SX7v x ) + -^-{uj<j)S<t> + 5F- HSF)S . (46) 

a * zF 

In the following we consider a situation in which the scalar-field dependent terms on the r.h.s. of Eq. (|45[) is neglected 
relative to the term Av x /a. In this case Eq. (|46p yields 

8 - k ~ -V • (S\7v x ) . (47) 

Later we shall confirm the validity of this approximation. 
From Eq. flSS} with Eqs. dM]), J35D and (JSHJ) we obtain 



/ F \ 1 

I 2Fy 2F 

iV F 3F . 
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dp + 4w4>8<j> + {2uj^4>' 2 + FjR - 2V^)84> + 38F + 3HSF -\&H 2 + —)8F 
3H + ^-(6F + 6HF + Aujtf) + -5. 



(48) 



where Nn and TVs correspond to the second-order terms on the r.h.s. of Eqs. (j29f and (I32)) . respectively. Following 
Refs. [3j, [y, |27I - |45j we employ the sub-horizon approximation under which the terms containing k, 8p, A 8F/a 2 and 



Aa/a 2 are picked up in Eq. (|48|) . Note that \F/HF\ <C 1 under the condition \<j)\ <C H. Apart from the term Aa/a 2 , 



the terms on the r.h.s. of Eq. 



are of the order of H 2 a or smaller. We then have 



2F«- -L (*,-*«■ 



- 2 [(vv x ) ■ (Vv x y] 



(49) 



Of course this approximation is justified when the second-order term on the r.h.s. of Eq. (|49p is larger than the 
first-order (field-dependent) terms on the l.h.s. of Eq. (|48p we have neglected. Later we shall derive conditions under 
which this approximation is valid. 

Let us estimate the field perturbation 5<p as well as SF. As we explained in the previous section, the field mass 
M defined in Eq. (j2"Tj) is much larger than H. Using the approximation in which the terms containing A/ 2 , A 
ASF/a 2 , Sp and SR are dominant contributions to Eqs. (j3"Tj) and (|55|) . we obtain 



A 

~n2 






~Av x --F,SR: 



0. 



A 1 F 

SF--6p--Av x 



l FSR 



0. 



(50) 



(51) 



This approximation is accurate as long as an oscillating mode of the field perturbation does not dominate over the 
matter-induced mode |g,l22t- Note that we have neglected second-order terms on the r.h.s. of Eqs. (j3"7| and (|38[) . 
Since the field is nearly frozen at the instantaneous minimum given in Eq. (|18|). the dominant second-order term 
corresponds to V^^Scj) 2 . This term gives rise to only a tiny correction to the growth rate of perturbations. Moreover 
it can be neglected relative to the second-order term on the r.h.s. of Eq. ([49]) . See Appendix for the detailed estimation 
of such a second-order term. 

On combining Eqs. (|50|) and (J5T|) . we find 



M 2 
~F~ 



A 



5F 



2Q 2 Sp - 



F 



-Av^ 



(52) 



Note that 5p — p m S ~ 3FH 2 S during the matter era. At the first-order we also have the following relation from 
Eq. g7]): 



a x 



5 = cHS , c = D/HD , 



(53) 



where D(t) is the time-dependent part of S. Since D(t) is typically proportional to t n with n of the order of unity 
[13, it follows that c = 0(1). Hence we get 



(F/a)Av, 



2Q 2 Sp 
As long as the condition, 

1^1 « \QH\ 
is satisfied, we have that \2Q 2 Sp\ 3> \(F/a)Av x \ and 

'M 2 A 



QH 



F 



SF ~ 2Q 2 p m S . 



(54) 



(55) 



(56) 



In the previous section we showed that \<j>\ is much smaller than H for the potential @. Hence the condition ([55]) 
holds well for the values of \Q\ which are not very much smaller than 1. 

Equation (|56|) shows that the perturbation SF is sourced by the matter perturbation 6. Hence Eq. (|49|) can be 
written as 



2Hk - 4irp m G cS 5 



1 



[(Vv x ) ■ (vv x y 



(57) 



where 



G c ft8 = 



8nF 



8- 



1 ASF\ 



Pm 



a 2 J 



(58) 



We introduce an effective potential $ and a peculiar velocity u as follows: 



a 



2 



iTTpmGesS , (59) 

-Vv x . (60) 



If one defines an effective gravitational potential $ = <p + aHv x in Eq. ([30]) . it follows that A^/a 2 = 4irG c ffp m 5 at 
linear order, where G c ff is different from G Q g by the sign of AF/a 2 . Taking the time-derivative of Eq. (|47|) together 
with the use of Eq. (|57j) . we get 

<9 2 <S 9<J 1 „ , 1 d 2 , , ,, , N 

a? + 2H m = ^ ' (1 + * )v$ + ^a^ < uV ) ■ (61) 

This is our master equation that is used to compute the skewness of matter density perturbations in the next section. 

IV. SKEWNESS IN MODIFIED GRAVITY 

In this section we study the skewness of matter perturbations for the action ([2]) with the potential (|4]) . The skewness 
will be derived analytically in the matter-dominated epoch by using Eq. (foil) . 

A. First-order perturbations 

We write the solution to Eq. (|6"Tj) in the form 5 — <5 < - 1 - ) + 5^ + ■ ■ ■ , where the subscripts represent the orders of 
perturbations. The equation for the first-order perturbation S^ 1 ' is 

^ + 2I^-4W^=0. (62) 

We express the first-order perturbations 5^ and SF^ in the plane-wave form: #W = J Sj, (t)e~ lkx d 3 fc and 8F^ = 
J^ 1) (t)e- ifea5 d 3 fe. From Eq. (JSSJ) we obtain 

^"'"' Mvmv/ '"'"' (63) 

Then the temporal part of Eq. (|62|) satisfies 

6£\t) + 2HS { k 1 \t) - Anp m G^8^(t) = , (64) 



where 



,(i) _ I (k 2 /a 2 )(l + 2Q 2 )+M 2 /F 
8irF (k 2 /a 2 )+M 2 /F 



G% = ^ ' nl^C LZZ ' ■ (65) 



In the early stage of the matter era, the mass M is sufficiently heavy to satisfy the condition M / F 3> k /a . In 
this regime we have G e ^ ~ 1/8ttF ~ G, thus mimicking the evolution in General Relativity. At late times it happens 

that the perturbations enter the regime M 2 /F <C k 2 /a 2 . This case corresponds to G C J ~ (1 + 2Q 2 )/8nF, thus 
showing the deviation from General Relativity. The transition from the former regime to the latter regime occurs at 
a redshift Zk given by [27| : 



Zk 



]_\2(1 P ) 2 p pC 1 y^ 



a H Qj {l-pY- p {iF n^) 2 -pH 2 



1 p 

- 1 , (66) 



where ao and Hq are the present values. 

Using the derivative with respect to N = In (a), Eq. ([64]) can be written as 



6 * } " + \\- l WcS ) 5 * ] ' 12 ^ FQ mG^5^ = , (67) 
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where a prime represents the derivative in terms of N. As we explained in Sec. HH we are considering the case in 
which the matter era is realized by the point (M2) with |A| 3> \Q\ = 0(1). Since Sl m ~ 1 and w e g ~ in this case, 
we get the following solutions 



?(i 



'i 2 / 3 . 



for z » Zk, 



6 k ( i )«i^(V^i^-i) ; forz<<Zfe 



(68) 



In these asymptotic regimes the growth rates of first-order perturbations are independent of the wavenumber k. The 
growth rate is constrained to be s = S k /H8 k < 2 from observational data, which gives the bound \Q\ < 1 [271 ] - 

B. Conditions for the validity of approximations to reach Eq. (l(il[) 

In order to reach Eq. (|47p . we have employed the approximation that the field-dependent term on the r.h.s of 
Eq. (jlij)) is neglected relative to the term Av x /a. We have also neglected some of first-order terms in Eq. (jig)) relative 
to the second-order term 4j- [(Vw x ) • (Vi> x )' z ] . in Eq. (|49[) . We derive conditions under which these approximations 
are justified. 

We write the temporal part of the first-order perturbation v x as (v x ) k (t). Using Eqs. ([63]) and ([53]) together with 
the relation p m ~ 3FH 2 that holds during the matter era, we find 



HSF, 



(i) 



Q 2 H 2 



M 2 /F + k 2 /a 2 



HS 



(i) 



Q 2 H 2 



M 2 /F + k 2 /a 2 



V*)L 1} 



(69) 



As we showed in Sec. [II] the condition, M 2 jF ^S> H 2 , holds for the potential (HJ. Moreover we are considering sub- 
horizon modes deep inside the horizon, i.e., k 2 ^ a 2 H 2 . This leads to the relation \HSF k / F\ <C \(Av x ) k /a\ in 
Eq. (|69[) . Since the field (f> is nearly frozen at instantaneous minima of its effective potential, we have the relation 

\6F, 



frCOi < 



rWl 



m 



\(» 



\H5F^ '\ and hence \5Ffc ' /F\ <C \(Av x ) k /a\. The following inequality is also satisfied: 



to^ k 1] 



(1 - 6Q 2 )0 SF t 



(i) 



2Q 



< 



(Arf 



(70) 



where we used Eq. ([55]) . The above discussion shows that the field-dependent terms in Eq. ([15]) are neglected relative 

to the term Av x /a, thus ensuring the validity of the approximation, iv- 1 ' ~ Av x ' /a, for the modes deep inside the 
Hubble radius. 



The second-order term on the r.h.s. of Eq. (|4"9"]l is of the order of H 2 



([55]) . Meanwhile one of the first-order term, H 2 8F k / F, in Eq. 
former is larger than the latter provided that 



?(!) 



by employing the first-order solution 



has been already estimated in Eq. 



i4 1} i» 



Q 2 H 2 



M 2 /F + k 2 /a 2 



The 



(71) 



The r.h.s. of Eq. ([7T]) is much smaller than unity because of the condition {M 2 /F, k 2 /a 2 } ^> H 2 . One can show that, 
under the condition (j7lT) . other field-dependent first-order terms on the l.h.s. of Eq. ([48]) can be negligible relative to 
the term -^ [(Vw x ) • (Vu x ) : *l .. In summary, the master equation (|6ip we have approximately derived in the previous 
section is trustable under the conditions ([55]) and ([7T]) . 



C. Second-order perturbations and skewness 



The second-order perturbation 5^ satisfies 
d 2 6^ _. d5& 



dt 2 



2H d -^- - A, Pm G^ = 47rG « Pm (,(D) 2 + l_ d (?$f + J_ [ 



u iw u m 



(72) 



where G cfl 8^ = [S^ — ASF^ /(p rn a 2 )]/(87rF). When the growth rate of perturbations is dependent on k, the 
gravitational constant G C J is generally different from G ci f . In the following we study the regime of the massless limit 
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[M — > 0) in which the growth rate of the first-order perturbation is independent of k, i.e, S 1 -- 1 ^ = D{t)5\{x) with 
D(t) = t(\/25+48Q^-i)/6_ In this regime we h ave £(2) = (jW = (1 + 2Q 2 )/8nF, so we simply adopt the notation, G cS , 
instead of G ef l and G cS ■ . The General Relativistic case (M — * oo) is recovered by taking the limit Q — v 0. 



The first-order solution to u can be obtained by solving Eq. ([55]) , i.e., V • m^ 1 ) = — a5^\ It is given by 

a) aD f (x — x')8i(x') ,, , at) 



47T 

where A j2 ; is a spatial derivative of the quantity: 

A(x) 

This satisfies the relation A^ = —4ir5i(x). 
The last term on the r.h.s. of Eq. ([72]) yields 



„/|3 



d 3 x > = rrlA 

47T 



<y i( !E ') j3„/ 



as — x' 



7 d J a;' 



(73) 



(74) 



1 

~n2 



u m u m 



l r 



D 2 



f ) 4 1) + 24! ) ^ (1) + 4 1) ^ 1) 






16ir 



2 A^jA^j 



We write the solution of Eq. (f?2")) in the form [3l|] 



<5 (2) =4 2) +4 2) , 



where ft ^ and S^ ' satisfy 



d 2 si 2 \ asP 

9 Sb ' 2H^ - 4nG eS p m S™ = D 2 



D 



4irG eSPm 8^ = AnG oSPm D 2 8 2 + — ^j >6 hi 



Ot 2 



Of 



5\ - —6 lti A, 

Zir 



16tt 



2^,y^,'i 



(75) 
(76) 

(77) 
(78) 



Since $ 



(i) 



-G eS p m a 2 DA ti from Eq. (JSSJ, the r.h.s. of Eq. fl77|) is given by 4irG cS p m D 2 [5 2 - (l/47r)A ?i (5i 



s( 2 ) 



Writing the solution of ft as ft = E a (t)S a (x), we obtain the following equation for the temporal part: 

£„ + 2HE a - 4-KG eS p m E a = 4irG cS p m D 2 , 



(79) 



where the spatial part is given by 5 a (x) = S 2 — (I/AttjAjSij. Expressing the solution of Eq. (|T8"|) in the form 
6 { b 2) = E b (t)S b (x), we get 

E b + 2HE b - 4irG cSPm E b = D 2 , (80) 



and ft(x) =S 2 - (l/2ir)A ti 6 h i + (l/16vr 2 )A^A i4i . 

We then find the following solution for second-order perturbations: 



6^(t,x) = E a (t) 



Air 



E b {t) 



5\ - — A ti 6 lti 

Z7T 



Wir 2A ' ijA ' ij 



D 2 + E a 



D 2 
§1 - i-A,i<5 M 

47T 



D2 - E ^ A 



(81) 



In the second line we employed the fact that E a and 25& are related each other via the relation E a + 2E b = D 2 . 

We assume that the initial distribution of perturbations is Gaussian so that it is described by an auto-correlation 
function £(x) satisfying 

(*(»)) =0, (6(xi)d(x 2 )) = £(\xi - x 2 \) , (8(xi)S(x 2 )d{x 3 ))=0, 

(S(x 1 )S(x 2 )5(x 3 )5(x i )) = a\xi - x 2 |)£(|x 3 - x 4 |) + a\xi - x 3 |)£(|x 2 - x 4 |) + £(|xi - x 4 |)£(|x 2 - x 3 |) . (82) 
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Since ((S^) 3 ) = 0, the quantity (<5 3 ) is given by (S 3 ) = 3((<5W) 2 (5( 2) ) to the lowest order. We then have 

<<5 3 ) = \d 2 {D 2 + E a ){6$) - ^D*(5?5x, 4 A, 4 ) + ^D 2 (D 2 E a )(6 2 A^A, tJ ) . (83) 



Since the each ensemble average in Eq. (|83[) satisfies the relations (5f) = 3£(0) 2 , (Sfdi^A^) — 47r£(0) 2 and 
■A ■■) - 202C 



(SfA tij A tij ) = ^f-£(0) 2 |3fl|, we obtain the skewness 



53 "W = 4 + 2 §' (84) 

where we used (S 2 ) = £(0) 2 . Hence the skewness is determined by the second-order growth rate E a relative to the 
squared of the first-order growth rate D. 
Equation ([75)1 can be written as 

K +{\~ I w cb) K - 12irFG cB n m E a = \2vFG cS n m D 2 . (85) 

Recall that G c s = (1 + 2Q 2 )/8irF in the limit M — > 0. During the matter era realized by the point (M2) we have 
f2 m ~ and w c g ~ 0, in which case Eq. ([85| reduces l 

K + \K - \ (1 + 2Q 2 )E a = |(1 + 2Q 2 )i? 2 . (86) 



Using the first-order solution D — e4(v 25 + 48< 3 2 ^ N y we get the following special solution for Eq. (|8"6"| : 



19 + 36Q 2 - v/25 + 48Q 2 



Hence the skewness in the regime of the massless limit is given by 



s = 4[22 + 42Q 2 - ^25 + 48Q 2 ] 
19 + 36Q 2 - y/2b + 48Q 2 

The General Relativistic case is recovered by taking the limit Q — > 0: 

5 3 = 34/7 . (89) 

This agrees with the skewness in the Einstein-de Sitter Universe [301 (pressureless matter without cosmological con- 
stant). 

In Fig. Q] we plot the analytic value (|88j) as a function of |Q|. The skewness shows some difference compared to the 
Einstein-de Sitter value 34/7 for \Q\ > 0.1. When \Q\ — 1 we have S 3 = 4.775, which is different from the value 34/7 
only by 1.7 %. For the potential ([J| the first-order perturbation S k evolves from the regime M 2 /F ^> k 2 /a 2 to the 
regime M 2 jF <g; k 2 /a 2 for the modes relevant to large-scale structure. Hence the skewness tends to evolves from the 
value 34/7 to the asymptotic value given in Eq. ([55]) . The estimation (|55|) has been derived by neglecting the transient 
phase around the redshift Zk- Since this transition occurs quickly for the models that satisfy local gravity constraints 
(p > 0.7) [8|, [2(J, [27|, it is unlikely that the skewness is altered significantly by the presence of this transient phase. 

The estimation (|88| does not take into account the evolution in the late-time accelerated epoch. In the ACDM 
model the numerical analysis shows that the skewness increases a bit during the accelerated phase from the value 
34/7 (= 4.857) to the present value 4.865 (at il m = 0.28). This corresponds to the growth of 0.16 % only. We have 
checked that this situation does not change much even in the presence of the coupling Q. Hence the difference of 
the present values of the skewness from that in the ACDM model is only less than a few percent. This shows that 
the skewness provides a robust prediction for the picture of gravitational instability from Gaussian initial conditions, 
including scalar-tensor models with large couplings (|Q| < 1). 



1 Note that the skewness was calculated in Ref. |39( when the matter era is realized by the point (Ml). As we already mentioned, the 
point (Ml) can not be used for the matter era when the coupling |Q| is of the order of unity. 
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FIG. 1: The analytic estimation (|88p of the skewness during the matter-dominated epoch in the regime of the massless limit 
(M — » 0). With the increase of \Q\ the skewness gets smaller compared to the value 53 = 34/7 in the Einstein-de Sitter 
Universe. However the difference from the Einstein-de Sitter case is only less than 1.7 % for \Q\ < 1. 



CONCLUSIONS 



In this paper we have studied the evolution of second-order matter density perturbations in a class of modified gravity 
models that satisfy local gravity constraints. We have considered the scalar-tensor action @, which is equivalent to 
Brans-Dicke action (JXJ) with the correspondence 1/(2Q 2 ) = 3 + 2cjbd- In the presence of a field potential it is possible 
to satisfy local gravity constraints (LGC) even when \Q\ is of the order of unity. In fact the potential ((4]) is designed 
to have a large mass in the region of high density for the consistency with LGC. This covers the models proposed by 
Hu and Sawicki [16| and Starobinsky [8[ in the context of /(i?) gravity (Q = — 1/V6). 

Starting from second-order relativistic equations of cosmological perturbations, we have derived the equation (frJTj) of 
matter density fluctuations approximately. In so doing we employed the approximation that first-order perturbations 
in the scalar field (f> is neglected relative to second-order matter and velocity perturbations. This is valid under the 
conditions (|55|) and (fTTj) , both of which can be naturally satisfied for the values of Q we are interested in (0.1 < \Q\ < 1). 
Compared to the ACDM model, the effective gravitational constant G c g is subject to change at the late epoch of the 
matter era. This leads to the larger growth rate of first-order matter perturbations (S k oc tW 25 + 48 Q - 1 )/ 6 ) compared 

to the standard case (S k ex £ 2 / 3 ). 

The skewness of matter distributions is determined by the second-order growth factor E a relative to the squared of 
the first-order growth factor D. In the "scalar-tensor regime" where the effective gravitational constant is given by 
G fi ~ (1 + 2Q 2 )/8nF, we have derived the analytic expression ([55)1 of the skewness in the matter-dominated epoch. 
In the "General Relativistic regime" where G c ff — l/87rF ~ G, we have reproduced the standard value S3 = 34/7 
in the Einstein-de Sitter Universe. In modified gravity models with \Q\ < 1, the analytic value ([55)) of the skewness 
in the asymptotic regime of the matter era is different from the value 34/7 only less than a few percent. Even if we 
take into account the evolution of perturbations during the accelerated phase, the difference of the skewness relative 
to the ACDM model remains to be small. The above result comes from the fact that the ratio of the second-order 
growth rate relative to the first-order one has a weak dependence on the coupling Q. 

When \Q\ = 0(1) the growth rate of first-order matter perturbations is significantly different from that in the 
ACDM model. This gives rise to large modifications to the matter power spectrum as well as to the convergence 
spectrum in weak lensing, while the skewness is hardly distinguishable from that in ACDM model. This fact can 
be useful to discriminate large coupling scalar-tensor models among many other dark energy models from future 
high-precision observations. 
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Appendix 

In this Appendix, we estimate the order of second-order terms that we have neglected in Eqs. ([50]) and l|51jl . The 
dominant contribution of such second-order terms comes from the third-derivative of the potential, i.e., V^^Stp 2 [44| . 
Compared to the field-mass dependent term V^Stf* on the l.h.s. of Eq. (|50p . we have {V^^Sif) 2 ) / {V^S^) w —S(f>/(j) 
for the potential (J4j) under the condition \Q<j>\ <C 1. 

Since the field stays around the instantaneous minimum given in Eq. (|18[) . the field 4> can be estimated as 

<f>~ 3(1 -p)Q^L, (90) 

where we used Eq. (|2"Tj) . In deriving this, we have also employed the approximate relation p m m 3H 2 . Note that the 
order of p m is not different from 3H 2 even at the present epoch. Meanwhile, from Eq. (|63[) . the first-order perturbation 
5(f> k in the Fourier space during the matter era is given by 

6< Pk -~M2,l.2/ n 2 d k ■ I 91 ) 



Hence we obtain the ratio 



M 2 + k 2 /a 2 



£' i M2 ,(d 



C , (92) 



<j> l-pM 2 + k 2 /a 2 k 

which shows that \5<j>£ /$ < 1 for 5^ < 1. 

The presence of the second-order term V^^84> 2 gives rise to a correction of the order M 2 8(j) k ' /<f> to the mass 
squared M 2 in Eq. (I55|) . In two asymptotic regimes (i) M 2 3> k 2 /a 2 and (ii) M 2 <C k 2 /a 2 , this appears only as 
next-order corrections to the small expansion parameters (k 2 /a 2 )/M 2 [regime (i)] and M 2 /(k 2 /a 2 ) [regime (ii)]. 

In the regime M 2 3> k 2 /a 2 the correction from the term V^^Scj) 2 to the effective gravitational constant G e J is 

estimated as 5G ci { rs Q 2 (k 2 / a 2 M 2 )8<j) k /<f>. This gives the correction to the third term on the l.h.s. of Eq. (fST| in 
the Fourier space: 



47rp m SGi^ 4 X) 



Q 2 fc 2 /« 2 ^ 2f (D 

1 - p M 2 



H 2 8 k 1] , (93) 



which is much smaller than the second-order term on the r.h.s. of Eq. (|57|) that is of the order of H 2 5 k 
In the regime M 2 < k 2 /a 2 we have 6G^ w Q 2 M 2 (a 2 /k 2 ) 64 > k 1} /<j> and hence 



.(i) ^(i) 



^ 2 ^. A (D 2 



4^^-^ -r^^pT^) ^r . m) 



which is again much smaller than the r.h.s. of Eq. ()57j) . 

The above estimation shows that neglecting second-order terms in Eqs. (|50|) and (|5ip is justified. 
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